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Abstract 

Multi-channel  Rayleigh  wave  data  were  acquired  across  aim  diameter  spillway 
tunnel  along  three  parallel  traverses  with  surface  to  tunnel  separations  of  0.90  m,  2.15  m 
and  3.13  m  depth,  respectively.  These  surface  wave  data  were  acquired  by  placing  a  24- 
channel  geophone  array  perpendicular  to  the  center-line  of  the  spillway  tunnel  and 
incrementally  moving  the  array  across  the  tunnel.  The  near  source-receiver  offset  was  6 
m;  the  4.5  Hz  geophones  were  spaced  at  0.5  m.  The  tunnel  locations  were  identified 
visually  on  velocity-filtered  common-shot  gathered  field  records.  Tunnel  locations  were 
also  identified  by  analyzing  common  shot-gathered  records  using  two  newly-developed 
automated  interpretation  programs:  Spiking  Filter  Analysis  and  Attenuation  Analysis  of 
Rayleigh  Waves  (AARW).  Electrical  resistivity  data  was  acquired  along  each  traverse  for 
comparison  purposes. 


Introduction 

The  engineering  geophysics  community  has  recently  focused  on  the  use  of 
Rayleigh  (surface)  wave  methods,  such  as  Multichannel  Analysis  of  Surface  Waves 
(MASW)  and  Refraction  Micrometer  (ReMi),  to  detect  manmade  tunnels  in  the  earth’s 
shallow  subsurface  (Miller  et  al.,  2006).  Successful  tunnel-detection  applications  of  these 
methods  have  been  reported.  However,  further  experimental  and  analytical  investigations 
are  required  to  comprehend  all  significant  aspects  of  the  observed  surface  wave  data. 
This  case  study  reports  on  three  alternate  surface  wave  methods  that  were  used  to  locate  a 
1  m  diameter  tunnel:  visually- identified  diffracted/reflected  surface-wave  energy;  the 
Attenuation  Analysis  of  Rayleigh  Wave  (AARW);  and  Spiking  Filter  Analysis. 

The  study  site  is  located  at  Ber  Juan  Park  in  the  City  of  Rolla,  Missouri  about  172 
kilometers  (107  miles)  southwest  of  St.  Louis  off  of  U.S.  Interstate  Highway  44. 
According  to  design  specifications,  the  tunnel  is  a  reinforced  concrete  lined  spillway 
pipe,  1.07  m  in  diameter  with  a  1.5  percent  gradient  downstream.  The  reinforced 
concrete  liner  is  0.06  m  thick  with  an  internal  diameter  of  0.94  m.  The  spillway  was 
sealed  using  a  compacted  high  plastic  clay  backfill  around  the  pipe  to  0.25  m.  The 
earthen  dam  containing  the  spillway  pipe  is  composed  of  compacted  local  clayey  fill  with 
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6  to  12  inch  stone.  The  spillway  pipe  did  not  contain  water  at  the  time  of  the  geophysical 
survey. 


The  primary  objective  of  this  case  study  was  to  determine  the  utility  and 
reliability  of  the  automated  Spiking  Filter  and  AARW  methods  for  shallow  tunnel 
detection.  The  tunnel  locations  output  by  the  automated  interpretations  were  compared  to 
the  spillway  blueprints  to  assess  accuracy.  Electrical  resistivity  profiles  were  acquired 
along  the  surface-wave  traverses  for  comparison  purposes. 


Figure  1:  (a)  4.5  Hz  geophone  array  over  spillway  tunnel  at  a  Ber  Juan  Park,  Rolla,  MO.: 
Graduate  student  Wamweya  uses  a  sledge  hammer  source  to  generate  surface  wave 
energy,  (b)  Design  blueprints  and  on-site  measurements  of  the  Ber  Juan  spillway  pipe  are 
compared  to  the  geophysical  interpretations  to  assess  the  reliability  of  automated 
interpretation  methods. 


Resistivity 

An  AGI  SuperSting  automated  resistivity  system  was  used  to  image  the 
subsurface  across  the  spillway  tunnel  along  the  surface  wave  traverses  (Figure  2).  A 
dipole-dipole  array  was  employed.  Electrodes  were  spaced  at  0.5  m  (for  the  0.90  m  and 
2.15  m  depths  tests)  and  at  1  m  (for  the  3.13  m  depth  test).  The  air-filled  tunnel  is 
represented  as  a  zone  of  anomalously  high  resistivity  on  the  three  profiles  (Figure  2).  The 
resistivity  profiles,  as  a  whole,  exhibit  relative  homogeneity  below  the  thin  layer  of  dry 
capping  soil  Resistivities  in  the  range  of  15  to  20  Qm  are  consistent  with  clayey 
compacted  fill  (Reynolds,  2000;  p.  442).  The  interpretation  of  the  resistivity  data 
supports  that  the  subsurface  geology  at  Ber  Juan  is  effectively  uniform  (except  for  the 
tunnel)  and  therefore  an  ideal  field  laboratory  for  this  study. 


Figure  2:  Scaled  electrical  resistivity  profiles  acquired  across  the  Ber  Juan  Spillway 
tunnel  at  tunnel  depths  of  (a)  0.90  m,  (b)  2.15  m  and  (c)  3.13  m,  respectively.  Horizontal 
and  vertical  scales  are  in  meters. 


Surface  Wave  Analysis 


Data  Acquisition 

The  Rayleigh  (surface)  wave  data  were  acquired  using  a  Seistronix  RAS-24 
seismograph  and  24-channel  streamer  consisting  of  4.5  Hz  geophones  spaced  at  0.5  m.  A 
10  kg  sledge  hammer  source  with  a  near-offset  of  6  m  was  employed.  Surface  wave  data 
were  acquired  along  three  parallel  traverses  oriented  perpendicular  to  the  axis  of  the 
tunnel.  The  tunnel  was  crossed  at  depths  of  0.90,  2.15  and  3.13  m,  respectively.  ReflexW 
seismic  software  was  used  for  pre-processing,  editing,  visualization  and  1D/2D  spectrum 
analysis.  The  average  Rayleigh-wave  phase  velocity  of  the  upper  3  m  of  shallow  soil  was 
—160  m/s. 


Common  Shot  Gathers 

The  location  of  the  tunnel  (immediately  beneath  trace  6)  can  be  estimated  from 
the  analyses  of  the  common-shot  gathered  field  data  (Figure  3).  More  specifically,  the 
shallowest  surface  waves  are  significantly  lower  amplitude  to  the  right  of  trace  6  than  to 
the  left  as  a  result  of  tunnel-related  attenuation.  Also,  backscattered  (reflected/diffracted) 
surface  wave  energy,  characterized  by  negative  apparent  velocities  (reverse  dip),  is 
identified  on  traces  1-5. 


Figure  3:  (a)  2D  Rayleigh  wave  profile  above  the  known  Ber  Juan  spillway  tunnel  at 
0.90  m  embedment  depth;  (b)  at  2.15  m  depth;  (c)  at  3.13  m  depth.  The  tunnel  is  located 
immediately  below  trace  6.  The  shallowest  surface  waves  are  significantly  lower 
amplitude  to  the  right  of  trace  6  than  to  the  left  as  a  result  of  tunnel-related  attenuation. 
Also,  backscattered  (reflected/diffracted)  surface  wave  energy,  characterized  by  negative 
apparent  velocities  (reverse  dip),  is  identified  on  traces  1-5. 


Velocity  Filter  to  Enhance  Back-Scattered  Energy 

Surface  wave  energy  is  reflected/diffracted  from  shallow  tunnels.  On  common 
gathered  data,  diffracted/reflected  surface  waves  are  characterized  by  hyperbolic  arrival 
times  (Figure  4).  To  accentuate  diffracted/reflected  surface-wave  energy,  an  F-K  filter 
was  applied  to  the  common  shot  gathered  data  (Figure  4).  Xia  et  al.  (2006)  present  a 
simple  method  to  detect  and  characterize  voids  directly  from  a  filtered  common  shot 
gather. 


Figure  4:  (a)  F-K  filtered  common  shot  gathered  data  (a)  tunnel  at  0.90  m  embedment 
depth;  (b)  at  2.15  m  depth;  (c)  at  3.13  m  depth.  The  locations  of  the  tunnel  are  coincident 
with  the  apex  of  the  respective  arrival-time  hyperbola. 

To  illustrate  Xia’s  approach,  consider  Figure  4c.  The  depth  to  the  top  of  the 
tunnel  can  be  estimated  using  the  two-way  travel  time  t0  (to  apex  of  the  hyperbolic;  0.065 
ms),  the  arrival  time  of  the  hyperbolic  event  at  trace  2  (0.080  ms),  and  the  average  phase 
velocity  of  the  shallow  soil  (160  m/s).  Although  this  approach  appears  to  be  very 
reasonable,  the  output  depth  estimates  are  not  always  as  reliable  as  one  might  hope. 
Using  the  above  input  values,  the  depth  to  the  top  of  the  tunnel  imaged  on  Figure  4c  was 


estimated  to  be  1.8  m,  whereas  it  is  known  to  be  3.13  m.  Giles  et  al.  (2005)  in  his 
numerical  modeling  results,  reports  that  cavities  with  circular  sections  generate  less 
diffraction  than  rectangular  sections. 

Spiking  Filter  Analysis  on  Common  Shot  Data 

The  spiking  fdter  analysis  used  in  this  study  is  described  as  an  algorithm  that 
compares  the  measure  of  energy  response  of  the  geophones  nearest  the  source  to  the 
remaining  geophones  in  the  array.  In  essence,  the  common  shot  gather  is  deconvolved  to 
the  first  geophones.  If  there  is  no  anomaly  within  the  region  between  the  excitation  pulse 
and  the  last  geophone,  the  spiking  filters  designed  on  the  first  few  geophones  will 
deconvolve  the  responses  from  all  the  geophones.  However,  if  an  anomaly  (such  as  a 
void)  does  exist,  the  resulting  reflections  and  diffractions  will  disturb  the  geophone 
responses  such  that  a  spiking  filter  designed  from  one  response  will  be  unable  to  properly 
deconvolve  the  others.  We  form  a  so-called  spiking  statistic  from  the  ratio  of  energies 
after  and  before  application  of  the  spiking  filters  for  each  shot.  A  high  spiking  statistic 
indicates  there  is  an  anomaly  present  in  the  shot  preventing  effective  deconvolution,  thus 
indicating  the  presence  of  a  void. 

Figure  5  below  shows  the  results  of  the  algorithm  applied  to  data  collected  where  the 
geophone  spacing  was  0.5  meters.  As  the  geophone  array  was  dragged  across  the  area 
oriented  perpendicular  to  the  tunnel,  shots  were  taken  at  0.5  meter  intervals.  The 
difference  in  energy  is  recorded.  Subsequent  common  shot  gathers  build  a  virtual 
common  offset  profile.  The  figure  depicts  excitation  over  the  tunnel  at  the  7th  shot,  the 
first  array  geophone  was  over  the  tunnel  on  the  13  th  shot  and  the  last  array  geophone  was 
over  the  tunnel  on  the  37th  shot. 


Figure  5:  ‘Spiking’  statistic  peaks  are  noted  just  after  the  excitation  passes  over  the 
tunnel  at  shot  7.  The  array  spacing  was  0.5  m  and  offset  3  m. 


Attenuation  Analysis  of  Rayleigh  Waves  (AARW) 

Numerical  studies  showed  that  a  rectangular  void  starts  to  vibrate  in  response  to 
the  Rayleigh  wave  excitation  and  due  to  the  void  vibration  energy  partitioning  occurs 
(Nasseri-Moghaddam  et.  al.  2006).  Part  of  the  incident  energy  is  reflected  in  the  form  of 
Rayleigh  wave,  another  part  is  converted  to  body  waves  and  spreads  into  the  medium. 
Another  part  of  energy  is  trapped  in  the  void  region  and  bounces  back  and  forth  between 
the  void  boundaries,  until  it  damps.  Therefore,  the  transferred  part  of  the  energy  is 
attenuated  and  has  smaller  amplitudes.  The  trapped  energy  is  associated  to  higher  modes 
of  Rayleigh  waves  and  excited  Lamb  waves.  The  effects  of  reflected  and  trapped  energies 
are  seen  as  a  region  in  the  vicinity  of  the  void  with  concentrated  energy,  both  in  time  and 
frequency  domains.  The  extent  of  this  region  depends  on  the  void  size,  and  the  frequency 
content  of  the  incident  energy.  Thus,  in  some  cases  it  is  possible  to  correlate  the  size  of 
the  model  to  the  extents  of  the  region  with  energy  concentration. 

A  technique  known  as  Attenuation  Analysis  of  Rayleigh  Waves  (AARW)  was 
proposed  to  determine  the  location  of  a  void,  and  estimate  its  embedment  depth  (Nasseri- 
Moghaddam  et  al.  2005).  This  technique  is  based  on  the  observed  damping  effect  of  the 
void  on  the  surface  responses.  In  summary,  the  method  suggests  to  calculate  the 
cumulative  energy  of  the  responses  at  the  location  of  each  receiver  over  the  reliable 
frequency  range  (energy-distance  parameter  -  NED)  and  plotting  it  versus  distance.  The 
location  of  the  peak  of  the  graph  corresponds  to  the  void  location.  To  estimate  the 
embedment  depth  of  the  void  the  logarithmic  decrement  of  the  frequency  data  obtained 
from  subsequent  receivers  are  calculated.  A  summation  of  the  obtained  value  is  carried 
out  over  the  distance  (CALD  parameter)  and  the  obtained  parameter  is  plotted  verses 
frequency.  The  plot  fluctuate  up  to  a  certain  frequency  (cutoff  frequency)  that  its 
wavelength  is  associated  to  the  embedment  depth  of  the  void.  Nasseri-Moghaddam 
(2006)  presents  the  details  of  the  method  along  its  application  to  the  data  obtained  from 
numerical  models. 

Figure  6  shows  the  contour  plot  of  the  collected  data  over  Ber  Juan  tunnel  in  time 
domain.  The  horizontal  axis  represents  time  and  vertical  axis  shows  the  distance  from 
source.  In  this  case  the  distance  between  the  source  and  the  first  receiver  is  1  m.  The 
dashed  lines  show  the  projected  boundaries  of  the  void  on  the  surface.  The  measured 
Rayleigh  wave  velocity  before  the  void  is  about  170  m/s.  It  is  seen  that  larger  wave 
velocities  are  measured  after  the  void  (slope  of  the  events  are  larger).  Previous  studies 
using  numerical  models  and  rectangular  voids  showed  similar  results  over  the  void.  The 
nature  of  this  behavior  is  under  further  study.  In  the  time  domain  data  clear  indications  of 
the  void  such  as  reflections  or  amplitude  increase  in  the  void  region  are  not  seen.  Figure  7 
shows  the  same  data  in  frequency  domain.  These  data  are  multiplied  by  a  gain  function 
to  eliminate  the  effect  of  geometrical  damping.  The  effect  of  void  is  seen  in  the  form  of 
energy  concentration  before  and  close  to  the  void.  Further,  a  region  with  small  frequency 
amplitudes  (damped  region)  is  observed  after  the  void.  Contrary  to  the  previous 
experiences  (Nasseri-Moghaddam,  2006)  the  energy  concentration  before  the  void  -  that 
is  associated  to  the  reflected  waves  -  is  more  visible  than  the  one  over  the  void.  The  latter 
differences  between  this  figure  and  the  ones  obtained  from  numerical  models  can  be 


associated  to  the  differences  between  the  shape  of  the  tunnel  and  the  shape  of  the  void  in 
the  investigated  numerical  models. 
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Figure  6:  Contour  plot  of  the  data  collected  over  Ber  Juan  spillway  tunnel  at  0.90  m 
depth  in  time  domain.  The  dashed  lines  show  the  near  boundary,  centerline  and  the  far 
boundary  of  the  tunnel. 
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Figure  7:  Contour  plot  of  the  data  collected  over  Ber  Juan  spillway  tunnel  at  0.90  m 
depth  in  frequency  domain.  The  dashed  lines  show  the  near  boundary,  centerline  and  the 
far  boundary  of  the  tunnel. 


To  locate  the  void  the  normalized  energy  distance  parameter  (NED)  parameter  is 
calculated  and  plotted  versus  distance  (Figure  8a).  The  near  and  far  boundaries  of  the 
tunnel  are  shown  as  vertical  dashed  lines.  A  peak  with  a  sharp  decrease  in  the  value  is 
seen  right  at  the  near  boundary  of  the  void,  which  is  an  indication  of  the  location  of  the 
void.  However  due  to  the  existence  of  high  energy  concentrations  before  the  void  (in 
frequency  domain)  the  maximum  peak  of  this  curve  occurs  before  the  void  location. 
Figure  8b  shows  the  cumulative  amplified  logarithmic  decrement  (CALD)  versus 
frequency.  A  sharp  change  in  the  slope  of  curve  is  observed  at  approximate  frequencies 
85  Hz  and  170  Hz,  with  wavelengths  equal  to  2  m  and  1  m.  These  wavelengths  are 
associated  to  the  depth  to  the  top  and  bottom  of  the  tunnel. 


Figure  8:  a)  The  normalized  energy  distance  parameter  (NED)  verses  distance  from 
source,  b)  The  cumulative  amplified  logarithmic  decrement  (CALD)  value  with 
frequency.  Changes  in  the  slope  of  the  curve  are  observed  at  f  =  85  Hz  and  170  Hz,  which 
their  wavelengths  correspond  to  the  depth  to  the  top  and  bottom  of  the  tunnel. 


Summary  and  Conclusions 

Common  shot  gathered  Rayleigh  wave  data  were  acquired.  These  data  were 
visually  analyzed  after  an  F-K  filter  was  applied  to  enhance  diffracted/reflected  energy. 
AARW  and  Spiking  Filter  methods  were  applied  to  the  common  shot  data  to  estimate 
tunnel  parameters.  Complementary  resistivity  confirmed  the  encompassing  soil  was 
relatively  uniform.  The  known  dimensions  of  the  site  were  available  for  comparison 
through  design  drawings  and  on-site  tape  measurements. 

The  analytical  results  were  compared  using  a  MATLAB  routine  to  plot  the 
geospatial  estimates  of  tunnel  locations  (Figure  9b).  The  redundancy  of  using  this  multi¬ 
method  Rayleigh  wave  approach  reduces  the  risk  of  inconclusive  or  erroneous  data 
interpretation. 

Analyses  of  the  common  shot  gathered  data  (before  and  after  the  application  of  an 
F-K  filter)  illustrates  that  tunnels  at  the  Ber  Juan  site  can  be  visually  detected  on  field 
data  at  diameter  to  depth  ratios  as  low  as  1  to  3.  The  AARW  and  Spiking  Filter  methods 


appear  to  have  significant  potential  in  terms  of  the  automated  detection  of  tunnels.  In  our 
opinion,  further  work  to  quantifiably  characterize  tunnels  in  the  Earth’s  shallow 
subsurface  using  a  multi-method  Rayleigh  wave  approach  is  warranted. 
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Figure  9:  (a)  Laptop  with  Seismic  Interpretation  Software:  Here  the  operator  is  setting 
the  parameters  (e.g.  length  of  the  line  and  location  of  the  unit)  before  beginning  the 
survey;  (b)  Mapping  Tool  (Screen  Capture)  that  plots  anomaly  based  on  AARW  and 
Spiking  Analysis  where  certainty  is  based  on  relative  algorithm  strength  value  of  the 
anomaly. 
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Wood  Institute  from  the  U.  S.  Army  Research  Laboratory. 
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